An acoustic signature of extreme failure on model granular materials

Unexpectedly, granular materials can fail, the structure even destroyed, spontaneously in simple isotropic compression with stick-slip-like frictional behaviour. This extreme behaviour is conceptually impossible for saturated two-phase assembly in classical granular physics. Furthermore, the triggering mechanisms of these laboratory events remain mysterious, as in natural earthquakes. Here, we report a new interpretation of these failures in under-explored isotropic compression using the time-frequency analysis of Cauchy continuous wavelet transform of acoustic emissions and multiphysics numerical simulations. Wavelet transformation techniques can give insights into the temporal evolution of the state of granular materials en route to failure and offer a plausible explanation of the distinctive hearing sound of the stick-slip phenomenon. We also extend the traditional statistical seismic Gutenberg–Richter power-law behaviour for hypothetical biggest earthquakes based on the mechanisms of stick-slip frictional instability, using very large artificial isotropic labquakes and the ultimate unpredictable liquefaction failure.

For this study, the granular materials were created from industrial soda lime glass beads "Sil-glass" commercialized by CVP, Linselles, France [S1]. The quasi-perfect spherical glass beads SLG 6-8 have a mean grain diameter, d = D 50 = 0.723 mm, a mean surface roughness Sa, and a arithmetic mean height, of only 36 nm.
Short cylindrical granular sample, H 0 = 70 mm in height and D 0 = 70 mm in diameter, is fabricated using a modified moist tamping and under compaction method [S2, S3], mixed with 2% of distilled water by weight, then placed in the triaxial mould in five layers, compacted up to prescribed thickness inside an open ended cylindrical latex membrane of 0.3 mm of thickness.
To ensure a full saturated state, we use CO 2 method [S4] under very low gradient of only 0.2 kPa, followed by flushing de-aired and distilled water and application of constant back-pressure U 0 of 200 kPa to obtain a high Skempton's coefficient B = ∆U/∆σ 0.95 [S5].
For the simple isotropic compression under stress-controlled mode, we applied slowly compressed air to a desired cell pressure σ r = σ 3 , the radial (minor principal) stress. The sample is not kinematically constrained in the vertical direction since the vertical loading ram is not physically connected to the sample top cap. Consequently, the sample top cap is free to move vertically when sheared quasi-statically either in stress (isotropic compression) or strain-controlled (triaxial compression) loadings.
Comparing to past studies [S6, S7], the new acoustic measurements consist of one tiny piezoelectric accelerometer with resolution of 350 µg within 10 kHz (PCB 607A11) to measure the instantaneous vertical acceleration of the sample top cap inside the triaxial chamber; plus one complementary non-contact portable laser vibrometer, Ometron OH-100-D, and one free-field prepolarized microphone, BSWA MPA408 with ICCP preamplifier MA408, to measure the local lateral surface vibration at mid-height of the sample latex membrane and the radial sound pressure at the sample bottom outside the triaxial chamber.
We measure the pore water pressure at the sample top and bottom with two set of complementary pressure sensors: one static pore-water pressure U (Sedeme, MD20) with an unknown resonant frequency, complemented by one dynamic piezoelectric pore pressure sensor (PCB S112A21) with a high resonant frequency (250 kHz) to verify the fast rise-time and the unexpected oscillating nature of U .
We use two dynamic data acquisition boards NI 4472B of National Instruments to measure concurrently and synchronously all static and dynamic sensors at minimum 10 kHz to follow the fast dynamics of the isotropic collapse and liquefaction in detail.
The triaxial apparatus is located in a separated room and isolated from external building vibration by passive isolation systems. Fig. 2 gives a typical temporal evolution of isotropic collapse U 2 (blue) and liquefaction U 3 (red) in isotropic drained compression. A semilogarithmic time scale is used to follow the temporal behavior over fiver orders of magnitude with focusing on the short time scale below 10 ms. The time scale is shifted to the beginning of the changes in acceleration data ± 0.1 ms, which is the current time resolution. The first peak in the acceleration data, occurring at 0.8 ms and up to -1.5 g for the liquefaction event, is representative of a sudden fall of the sample top cap. This first acceleration peak always precedes that of the pore pressure by 2.0 ms. After this time delay, the pore fluid pressure changes rapidly from the imposed constant back pressure U 0 . Three phases characterize this sudden change [S6]. First, a transient phase I in which the pore pressure rises instantaneously from the constant back pressure U 0 and fluctuates toward a stabilised value U stable within 200 ms. Then an intermediate phase II at constant U stable before a slow dissipation phase III toward U 0 . The duration of phase II, longer than one second, is a characteristic for the specific liquefaction case, as U 3 . To compare different stick-slip like motions, the excess pore pressure U − U 0 is normalised by the stabilised value (U − U 0 ) stable .

B. Isotropic collapse and liquefaction dynamics
In saturated granular mechanics, liquefaction is associated with a vanishing effective stress, σ = σ -U = 0 [S8]; it means the total stress σ is equal to the pore fluid pressure U . For the liquefaction U 3 in Fig.1 in linear scale, the vertical red arrow indicates a null effective isotropic stress σ = 0 during a lengthy evolution of the solid fraction. For the same liquefaction U 3 in Fig.2 in logarithmic time scale, the normalized excess pore pressure, (U − U 0 )/(σ − U 0 ) = ∆U/∆σ = U/σ (for U 0 = 0 kPa) = 1 (red solid circle); thus σ = 0 and the liquefaction failure associated with a very large axial strain of 19.2% and a destruction of the sample shape.
Furthermore, for the liquefaction failure U 3 in isotropic compression under constant stress rate in Fig.1, we observe an unexpected jump in solid fraction of 0.018 (vertical red arrow at constant stress); and in the temporal complementary Fig.2, an exceptionally large incremental axial strain of 19.2% in one single step in less than one second. This duration of the liquefaction event of one second (or 0.037% of the total experimental duration) can be qualified as spontaneous regarding the usual time scale of the isotropic compression of about 45 mins to increase the effective cell pressure from 20 to 212 kPa.
The spontaneity of the failure event can be seen in the changes of the vertical acceleration in less than 0.01 seconds in Fig.2 (0.00037% of the total experimental duration), or in the changes of the pore fluid pressure in less than 0.1 seconds (0.0037% of the total experimental duration). Normally, these short-lived changes should not happen in a quasi-static regime under low dimensionless inertial number I ≈ 10 −6 [S9] in drained loading condition (constant pore fluid pressure).
Liquefaction in isotropic compression occurs under specific conditions [S6]: (1) the duration of phase II longer than one second, (2) the normalised excess pore pressure equals unity or the effective stress approaching zero value, and (3) an initially loose granular assembly above a threshold void ratio. If the duration of phase II is less than one second, the instability event is termed as isotropic collapse, and the compression experiment can resume without the destruction of the granular structure.
The incremental axial strain ∆ε 1 and the incremental volumetric strain ∆ε v show a dynamic axial slip in less than one second; and the the short-lived evolution of the pore fluid pressure happens essentially in undrained condition (no pore fluid flow). Fig. S1 gives the typical temporal evolution of the vertical top cap acceleration of the second event U 2 in the normal linear scale of 200 seconds from arbitrary beginning time, complementary to the logarithmic scale in Fig. ??. It shows a sharp spike of acceleration up to 1.5 g, the gravity acceleration g = 9.81 m/s 2 , at t = 11.571 seconds with practically no precursors (foreshocks) and very small aftershocks, comparing to real earthquakes. The top inset shows the close-up view of the main short-lived event U 2 of large magnitude near 1.5 g, lasting less than 5 ms; and the bottom inset the quasi-static regime of isotropic compression below 0.002 g under low stress rate ( 0.1 kPa/s) in a segment of 20 seconds between two main events U 2 and U 3 . The typical background noise level of acceleration signals is of about 29 dB. The cell pressure is constant at 295 kPa and the prescribed constant back pressure U 0 at 200 kPa in this segment.
The cylindrical sample shape is destroyed in a rare isotropic liquefaction failure U 3 , as shown in the supplemental movie M1a, acquired with fast camera at 3 000 fps (frames per second) and viewed at 30 fps. This movie is the first ever recorded isotropic liquefaction on fast camera. Fig. S2 (top) shows a typical spectrogram using a continuous Cauchy wavelet of 100th order of a quiet segment of one second of the vertical top cap acceleration preceding the second isotropic collapse U 2 , showing the nearly continuous frequency as a large horizontal red and blue band centered above ∼ 4 kHz. The close-up view in Fig. S2 (bottom) reveals the discrete details in a segment of 0.02 seconds (from 0.10 to 0.12 seconds) with a local short-lived power peak at 4.4 kHz and 0.115 seconds with a faint magnitude of only 6×10 −6 (au), about 4 orders of magnitude lower than that of typical isotropic collapse event. We can follow closely the temporal evolution of this local peak in the time-frequency analysis. Fig. S3 and table I give the Cauchy spectrogram and the vibrational characteristics by three complementary sensors in two perpendicular directions of the extreme failure event U 3 in a segment of 0.04 seconds, centred on the abrupt change of the vertical top cap acceleration G (top), the radial vibration with non-contact portable laser vibrometer V (middle) and the lateral sound pressure with free-field prepolarized microphone M (bottom). Each figure shows the temporal evolution of the time series in linear scale on top subfigure, its spectrum on the right subfigure and the spectrogram on the main panel. The vertical acceleration has the largest magnitude and the smallest duration, and two order of magnitude greater than that of lateral sound pressure. Fig. S4 indicates the time at peak amplitude and the very short total duration below 0.01 seconds for all isotropic events in a series of 17 experiments up to 500 kPa of confining pressure: collapse (a) and liquefaction (b), plus the liquefaction from the first stick-slip in triaxial compression (solid triangle) and Fig. S5 the associated dominant frequencies, especially the first mode around 300 Hz, function of stabilized excess pore pressure.

D. Multiphysics numerical simulations
We use Comsol Multiphysics [S10] with acoustic and geomechanic modules to explore the observed vibrational modes on a short full cylinder of homogeneous and simple plastic material consisting of isotropic linear elastic behaviour (Young modulus E, Poisson modulus ν and density ρ) with Mohr-Coulomb plastic criterion (cohesion c and frictional angle ψ). For simplicity, the two end conditions are clamped (no displacement for top and bottom caps), and free lateral surfaces are allowed. The finite element mesh of the thick cylinder has 14 000 tetrahedral elements; and that of the thin hollow cylindrical latex membrane enclosing the granular cylinder 15 000 tetrahedral elements.
With fixed ρ = 1.428 kg/m 3 , ν = 0.25 and cohesionless granular soil c = 0 kPa, ψ = 31 o , the first known measured modal frequency of 253 Hz for the liquefaction event U 3 at 212 kPa of confining stress is fitted using simple trial and error procedures to identify the Young modulus E = 11.74 MPa. It seems that the granular cylinder and the membrane are decoupled in numerical experiments since further simulations with a thin latex membrane enveloping a thick granular cylinder do not reveal the lower modes below 20 Hz.
In Fig. S6 on double logarithmic axis, the frequency of the vibrational modes and the estimated bulk modulus has a perfect power law relationship with an exponent of 0.5.

E. Acoustic energy
We use the temporal series of vertical top cap acceleration, as in Fig. 2, to determine the beginning of each instability event (isotropic collapse or liquefaction). The event signal was detrended, smoothed with Savitzky-Golay filter of 7th order to ease the event peak and duration detection. The local peaks are detected as the local maxima larger than the two neighboring data using the derivative based finding algorithm in "findpeaks" function of Matlab [S11]. Then the beginning t b(i) and the end t e(i) of each event are indicated by a sudden change of the first derivative of the smoothed temporal signal indicating the stabilisation backward and forward below a noise threshold from the local peak. The threshold for detecting the vertical acceleration is set to 0.005 g, larger than the actual noise level in Fig. S1. For isotropic liquefaction failure, the peak signal can be as large as 16 g.
The acoustic energy E a associated with the temporal series of the vertical top cap acceleration G(t) is estimated by the squared maximal value of G(t) with arbitrary unit. Regardless of the isotropic confining stress from 50 to 500 kPa, and of the uncontrolled triggering stress of the liquefaction failure in this stress range, the acoustic energy and the isotropic triggering stress, even uncontrolled, for the final liquefaction failure follow a highly correlated power-law (R 2 = 0.922) in Fig. S7 on double logarithmic axis.
The acoustic energy E ai can also be estimated by the integral of G(t) 2 over the event duration from the beginning t b(i) to the end t e(i) of the stick-slip event (equation S1). In our case, E a is proportional to E ai .
The probability distribution of E a is estimated on logarithmic bins with 10 bins par decade between the minimum and maximum of measured energy; and the frequency distribution is normalised by the bin width to give the probability density function (PDF). The power exponent is obtained using a nonlinear regression model "fitnlm" of Matlab based on iterative least squares Levenberg-Marquardt method.
F. Alternative analysis for the probability density function of acoustic energy The probability density function of the acoustic energy measured by the vertical top cap acceleration E a can be evaluated in a double power-law with a crossover energy E c ≈ 10 −5 (a.u.) in Fig. S8. Below E c , the energy distribution of small events yields a power exponent β 1 = 1.11 ± 0.04 (mean ± deviation, R 2 = 0.779) over only two decades, and another power exponent β 2 = 1.08 ± 0.01, R 2 = 0.951 beyond E c for nearly six decades. The smaller second power exponent β 2 of the double power-law, 1.08 instead of 1.21, seems to better suit the liquefaction data of large acoustic energy than the simple power-law. However, the relatively low correlation coefficient for β 1 makes this alternate analysis less preferable than the single power-law in the main text. Consequently, additional experimental works are needed to sort out the best analysis.
It seems that the lack of low energy cutoff in Fig. 5 is not a consequence of our experimental setup. The sensitivity of our accelerometer can measure very light energetic events, down to 10 −8 . Fig. S9 shows, as an example, the presence of a low energy cutoff E min = 2*10 −6 (a.u.) for water-wet model granular assembly in stick-slip experiment of triaxial compression at 500 kPa of confining pressure. Stabilized excess pore pressure ∆U stable (kPa) S5. Characteristic frequencies of collapse (a) and liquefaction (b) identified by CCWT from vertical top cap acceleration G, radial laser vibrometer V and lateral microphone M. SSLiq (StatLiq) means liquefaction from drained (undrained) stick-slip experiment.   . Double power-law relationship from:acoustic energy with vertical acceleration of the sample top cap : β = 1.11 ± 0.04, R 2 = 0.779 below Ec = 1*10 −5 (a.u.) and β = 1.08 ± 0.01, R 2 = 0.951 beyond Ec. The symbol ± stands for 95% confidence interval.